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.
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).
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 existsif (!dir.exists("soilgrid")) {dir.create("soilgrid")}# Create pathslfile <-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 layerssoilgrids_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.