TITLE: Querying the SoilGrids REST API
DATE: 2020-11-08
AUTHOR: John L. Godlee
====================================================================


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.

 [SEOSAW database]:

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.

 [SoilGrids]: https://www.isric.org/
 [REST API]: https://rest.soilgrids.org/soilgrids/v2.0/docs

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/laye
rs")
     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
 .         ...        ...      ...     ...