Querying the SoilGrids REST API

2020-11-08

This week I have been working on the SEOSAW database , which holds tree measurement data from around 7000 survey plots in southern Africa. The database is always difficult to deal with, mostly because the component datasets were all collected with different initial purposes, meaning there is a significant amount of data bashing to get each one into the SEOSAW format. I chose to use an R package to hold functions and workflows for cleaning the datasets and compiling the main database. In previous versions of the package I had a number of reduced size spatial data layers held as data objects in the package, which I could then query for information at each plot such as estimated herbivore biomass, elevation, and administrative region. I’ve since removed all of these functions, because they were too temperamental and the trade-off in spatial resolution that came from storing reduced size versions of each data layer in the package was too great. Now the data layers are queried at full resolution in a separate process right at the end of dataset compilation.

One of the functions which I was actually quite proud of was one that queried the SoilGrids REST API to get soil information. The function also provided functionality for averaging over soil depths. The SoilGrids REST API is still experimental, and since I first wrote the function they upgraded to v2.0. This broke the function and so it didn’t work for a long time.

I re-wrote the function yesterday and thought I would post it here, as it probably won’t get any use elsewhere. I also wrote a companion function which sends an empty query to get back all possible combinations of soil attributes, depths, and values (mean, Q0.05 etc.).

#' Query SoilGrids via the REST API 
#' 
#' See https://rest.soilgrids.org/soilgrids/v2.0/docs for more information
#' 
#' @param x dataframe of plot level data
#' @param plot_id column name string of plot IDs
#' @param longitude_of_centre column name string of plot longitude
#' @param latitude_of_centre column name string of plot latitude
#' @param attrib vector of soil attributes to extract
#' @param depth vector of soil depths over which to extract the attributes, e.g. 0-5, 30-60
#' @param average optional list of vectors, each of length 2, describing ranges 
#'     of depths over which to average (mean) the values of attrib
#' @param value vector of soil values to extract for each attribute, e.g. mean, Q0.05 
#' 
#' @details See \code{soilQueryOptions()} for all possibly combinations of 
#'     attrib, depth, and value
#'
#' @return dataframe of soil values by depth
#' 
#' @importFrom httr GET content
#' @export
#' 
soilQuery <- function(x, plot_id = "plot_id", 
  longitude_of_centre = "longitude_of_centre", 
  latitude_of_centre = "latitude_of_centre", 
  attrib = c("cec", "cfvo", "clay", "nitrogen", "ocd", "ocs", 
    "phh20", "sand", "silt", "soc"),
  depth = c("0-5", "5-15", "15-30", "30-60", "60-100", "100-200"), 
  average = list(c(0,30)), 
  value = "mean") {

  default_attrib <- c("bdod", "cec", "cfvo", "clay", "nitrogen", "ocd", "ocs", 
    "phh20", "sand", "silt", "soc")
  default_depth <- c("0-5", "5-15", "15-30", "30-60", "60-100", "100-200")
  depth_values <- unique(unlist(strsplit(default_depth, "-")))
  default_value <- c("Q0.05", "Q0.5", "Q0.95", "mean", "uncertainty")

  # Check input
  if (any(!attrib %in% default_attrib)) {
    stop("Illegal soil attribute: ",
      paste(attrib[!attrib %in% default_attrib], collapse = ", "),
      "\n\tAllowed values: ", 
        paste(default_attrib, collapse = ", "),
      "\n\tRun `soilQueryOptions()` to find valid attributes")
  }
  if (any(!depth %in% default_depth)) {
    stop("Illegal depth: ",
      paste(depth[!depth %in% default_depth], collapse = ", "),
      "\n\tallowed values: ", 
      paste(default_depth, collapse = ", "),
      "\n\tRun `soilQueryOptions()` to find valid depths")
  }
  if (!is.null(average)) {
    if (any(!unlist(average) %in% depth_values)) {
    stop("Illegal averaging depth: ",
        paste(average[!unlist(average) %in% depth_values], collapse = ", "),
        "\n\tallowed values: ", 
        paste(depth_values, collapse = ", "))
    }
    if (any(unlist(average) > max(depth_values) | unlist(average) < min(depth_values))) {
      stop("average range values outside sampled depths: ",
        paste(unlist(average)[unlist(average) > max(depth_values) | unlist(average) < min(depth_values)], 
          collapse = ", "))
    }
    if (any(sapply(average, function(y) {
          length(y) != 2
        }))) {
      stop("All average ranges must be vectors of length 2")
    }
  }
  if (any(!value %in% default_value)) {
    stop("Illegal value: ",
        paste(value[!value %in% default_value], collapse = ", "),
        "\n\tallowed values: ", 
        paste(default_value, collapse = ", "),
      "\n\tRun `soilQueryOptions()` to find valid values")
  }

  # Construct query
  attrib_string <- paste0(paste0("&property=", attrib), collapse = "")
  depth_string <- paste0(paste0("&depth=", depth, "cm"), collapse = "")
  value_string <- paste0(paste0("&value=", value), collapse = "")

  queries <- lapply(1:nrow(x), function(y) {
    call <- paste0("https://rest.soilgrids.org/soilgrids/v2.0/properties/query?", 
      "lon=", x[y, longitude_of_centre], 
      "&lat=", x[y, latitude_of_centre], 
      attrib_string, 
      depth_string,
      value_string)
  })
  
  # Run GET query
  query_get <- lapply(queries, httr::GET)

  # If any queries fail, end function
  if (any(unlist(lapply(query_get, `[[`, "status_code")) != 200)) {
    stop("Some queries failed")
  }

  # Flatten to list
  query_list <- lapply(query_get, httr::content, as = "parsed")

  # For each query, for each attrib, for each depth, for each value, extract values
  query_extract <- lapply(query_list, function(y) {
    lapply(y$properties$layers, function(z) {
      lapply(z$depths, function(v) {
        lapply(v$values, function(w) {
          w 
        })
      })
    })
  })

  # Make tidy dataframe of values
  extract_df <- data.frame(
    plot_id = rep(x[[plot_id]], each = length(attrib) * length(depth) * length(value)),
    attrib =  rep(rep(attrib, each = length(depth) * length(value)), times = length(x[[plot_id]])),
    depth =   rep(rep(rep(depth, each = length(value)), times = length(x[[plot_id]])), times = length(attrib)),
    value =   rep(rep(rep(value, times = length(x[[plot_id]])), times = length(attrib)), times = length(depth)),
    extract = unlist(query_extract))

  # Optionally average each value and attrib over depths
  if (!is.null(average)) {
    extract_split <- split(extract_df, 
      list(extract_df$plot_id, extract_df$attrib, extract_df$value))

    extract_df <- do.call(rbind, lapply(extract_split, function(y) {
      do.call(rbind, lapply(average, function(z) {
        y$min_depth <- sapply(strsplit(y$depth, "-"), `[`, 1)
        y$max_depth <- sapply(strsplit(y$depth, "-"), `[`, 2)
        out_df <- data.frame(plot_id = y[1,plot_id],
          attrib = y[1, "attrib"],
          depth = paste(y$min_depth[which(y$min_depth == z[1])], y$max_depth[which(y$max_depth == z[2])], sep = "-"),
          value = y[1, "value"],
          extract = mean(y$extract[which(y$min_depth == z[1]) : which(y$max_depth == z[2])],
          na.rm = TRUE))
      }))
    }))
    row.names(extract_df) <- NULL
  }

  return(extract_df)
}

The function to get all possible input options:

#' Query SoilGrids via the REST API 
#' 
#' See https://rest.soilgrids.org/soilgrids/v2.0/docs for more information
#'
#' @return Dataframe of all soil attributes, depths and values
#' 
#' @importFrom httr GET content
#' @export
#' 
soilQueryOptions <- function() {

  query_get <- httr::GET("https://rest.soilgrids.org/soilgrids/v2.0/properties/layers")
  query_list <- httr::content(query_get, as = "parsed")

  out <- do.call(rbind, lapply(query_list$layers, function(x) {
    attrib <- x$property
    attrib_df <- do.call(rbind, lapply(x$layer_structure, function(y) {
      depth <- y$range
      value <- unlist(y$values)
      data.frame(depth, value)
    }))
    cbind(attrib = attrib, attrib_df)
  }))
  
  return(out)
}

A test of soilQuery() for three locations in Africa:

x <- data.frame(plot_id = c("A", "B", "C", "D"), 
  longitude_of_centre = c(22.674 , 15.555, 32.122, -1.319738),
  latitude_of_centre = c(7.140, -14.736, -6.231, 53.065368))

attrib = c("cec", "clay", "nitrogen") 
depth = c("0-5", "5-15", "15-30", "60-100")
average = list(c(0,30), c(0,15), c(60,100))
value = c("mean", "Q0.5")

soilQuery(x, attrib = attrib, depth = depth, average = average, value = value)
plot_id attrib depth value extract
A cec 0-30 mean 102.00
A cec 0-15 mean 106.00
A cec 60-100 mean 100.00
B cec 0-30 mean 83.00
B cec 0-15 mean 84.50
B cec 60-100 mean 72.00
C cec 0-30 mean 53.00
C cec 0-15 mean 56.50
C cec 60-100 mean 44.00
D cec 0-30 mean 152.33
D cec 0-15 mean 173.00
D cec 60-100 mean 94.00
A clay 0-30 mean 187.00
A clay 0-15 mean 161.50
A clay 60-100 mean 335.00
B clay 0-30 mean 85.67
B clay 0-15 mean 73.50
B clay 60-100 mean 321.00
C clay 0-30 mean 108.00
C clay 0-15 mean 110.50
C clay 60-100 mean 265.00
D clay 0-30 mean 195.33
D clay 0-15 mean 211.00
D clay 60-100 mean 111.00
A nitrogen 0-30 mean 97.00
A nitrogen 0-15 mean 110.50
A nitrogen 60-100 mean 45.00
B nitrogen 0-30 mean 48.00
.