Accessing SoilGrids via {terra}

soil erosion

Author

Affiliation

Lomonosov Moscow State University

Published

March 5, 2022

After several QGIS and R packages updates, I cannot download SoilGrids with rgdal anymore. When I’m trying to run code from SoilGrids WebDAV tutorial I am receiving a following error:

ERROR 11: CURL error: SSL certificate problem: unable to get local issuer certificate

I never managed to fix this issue, so I decided to let it go since everyone is slowly quitting GDAL anyway. Here is my approach to downloading cropped and reprojected SoilGrids raster. There was a really great tutorial by Ivan Lizarazo on getting several SoilGrids layers using rgdal. My approach is mostly based on it.

1. Boundary layer

First of all, let’s load some boundary layer, i.e., our area of interest (AOI). I will use the default sf sample data for North Carolina counties. To reduce the size of AOI, let me select only the first county. The SoilGrids files are stored in Interrupted Goode Homolosine, so I have to reproject our AOI polygon to it.

Code

library(dplyr)
library(sf)

# Load sample data
nc <- st_read(
  system.file("shape/nc.shp",
    package = "sf"
  ),
  quiet = T
) %>%
  slice(1)

# Transform to IGH projection
igh <- "+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs"
nc_igh <- st_transform(nc, igh)

2. Download urls

Now Let’s just copy download urls from previous tutorials. And create a link to every other separate .vrt file, since we are gonna purrr::map them later. I’m interested only in mean topsoil characteristics (i.e. 0-30 cm) right now. So I will download sand, silt, and clay content and soil organic carbon content (soc).

Code

sg_url <- "/vsicurl/https://files.isric.org/soilgrids/latest/data/"
props <- c("sand", "silt", "clay", "soc")
layers <- c("0-5", "5-15", "15-30")

vrt <- paste0(
  props, "/",
  props, "_",
  rep(layers, 4),
  "cm_mean.vrt"
)

vrt[1]

[1] "sand/sand_0-5cm_mean.vrt"

Then, we need to create a list of paths to save. Let’s create a directory soilgrid where we are going to download our layers.

Code

# Optional
# Check if the directory exists
if (!dir.exists("soilgrid")) {
  dir.create("soilgrid")
}


# Create paths
lfile <- paste0(
  "soilgrid/",
  props, "_",
  rep(layers, 4),
  ".tif"
)

lfile[1]

[1] "soilgrid/sand_0-5.tif"

3. Download and preprocess function

My general idea is to crop the SoilGrid layer to the bounding box, reproject to my CRS (i.e. CRS of the nc layer), download, and then write as .tif. However, I want to do this for 12 rasters. Therefore, we need to write a function we are going to apply:

Code

library(terra)

# Function to download and transform soilgrid layers
soilgrids_download <- function(list_vrt, # download url
                               list_lfile, # destination path
                               shape_igh, # AOI shape in IGH proj
                               destproj) { # desired projection

  terra::rast(paste0(sg_url, list_vrt)) %>% # read raster
    terra::crop(ext(vect(shape_igh))) %>% # crop to bounding box
    terra::project(destproj) %>% # reproject
    terra::writeRaster(list_lfile,
      overwrite = T
    ) # Save
}

Before running it in the loop, let’s try it for the first layer.

Code

soilgrids_download(
  list_vrt = vrt[1],
  list_lfile = lfile[1],
  shape_igh = nc_igh,
  destproj = st_crs(nc)$proj4string
)

rast(lfile[1]) %>%
  plot()

It worked!

4. Map download

Next, with the help of purrr we can apply this function to all our links. Let’s measure elapsed time.

Code

library(purrr)
library(tictoc)

tic()
walk2(
  vrt,
  lfile,
  ~ soilgrids_download(.x, .y,
    shape_igh = nc_igh,
    destproj = st_crs(nc)$proj4string
  )
)
toc()

76 sec elapsed

Well, almost 2 mins. It is necessary to improve the timing somehow. This process can be run in parallel with furrr.

Code

library(furrr)

# Set Parallel
no_cores <- availableCores() - 1
plan(multisession,
  workers = no_cores
)

# Download!
tic()
future_walk2(
  vrt, lfile,
  ~ soilgrids_download(.x, .y,
    shape_igh = nc_igh,
    destproj = st_crs(nc)$proj4string
  )
)
toc()

44.31 sec elapsed

Code

# Exit parallel
plan(sequential)

Less than 1 minute. As Josh Starmer says, “BAM!”. Three times faster!

Code

list.files("soilgrid")

 [1] "clay_0-5.tif"   "clay_15-30.tif" "clay_5-15.tif"  "sand_0-5.tif"  
 [5] "sand_15-30.tif" "sand_5-15.tif"  "silt_0-5.tif"   "silt_15-30.tif"
 [9] "silt_5-15.tif"  "soc_0-5.tif"    "soc_15-30.tif"  "soc_5-15.tif"

Citation

BibTeX citation:

@online{tsyplenkov2022,
  author = {Tsyplenkov, Anatoly},
  title = {Accessing {SoilGrids} via \{Terra\}},
  pages = {undefined},
  date = {2022-03-05},
  url = {https://atsyplenkov.ru/posts/2022-03-05-soilgrids-terra/2022-03-05-soilgrids-terra.html},
  langid = {en}
}

For attribution, please cite this work as:

Tsyplenkov, Anatoly. 2022. “Accessing SoilGrids via {Terra}.” March 5, 2022. https://atsyplenkov.ru/posts/2022-03-05-soilgrids-terra/2022-03-05-soilgrids-terra.html.