Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Basic point-location REST requests to SoilGrids v2 API and wrangling of data into SoilProfileCollection object
# basic point-location REST requests to SoilGrids v2 API and wrangling of data into SoilProfileCollection object
# last update: 2020/08/18
library(httr)
library(jsonlite)
library(tidyverse)
# point id, latitude and longitude as inputs
your.points <- tribble(~id, ~lat, ~lon,
"A", 37.9, -120.3,
"B", 38.1, -121.5)
res <- lapply(split(your.points, f = your.points$id), function(yd) {
id <- yd$id
lat <- yd$lat
lon <- yd$lon
response <- httr::GET(sprintf("https://rest.soilgrids.org/soilgrids/v2.0/properties/query?lat=%s&lon=%s", lat, lon))
r.content <- httr::content(response, as = "text", encoding = "UTF-8")
res <- jsonlite::fromJSON(r.content)
extractSGLayerProperties <- function(x) {
out <- res$properties$layers[res$properties$layers$name == x,]$depths[[1]]
# fix names and labels for downstream
out <- out[,colnames(out)[grep("range", colnames(out), invert = TRUE)]]
out <- data.frame(label = gsub("cm", "", out$label), values = out$values)
colnames(out) <- gsub("\\.Q0\\.", "Q", colnames(out))
colnames(out) <- gsub("Q5", "Q50", colnames(out))
colnames(out) <- gsub("values", x, colnames(out))
return(out)
}
soc <- extractSGLayerProperties("soc")
bdod <- extractSGLayerProperties("bdod")
phh2o <- extractSGLayerProperties("phh2o")
clay <- extractSGLayerProperties("clay")
cec <- extractSGLayerProperties("cec")
# create new horizon data, merge in each property using depth label
hz.data <- tibble(id = id, lat = lat, lon = lon, label = soc[,"label"])
hz.data <- hz.data %>%
merge(soc, by = "label") %>%
merge(bdod, by = "label") %>%
merge(phh2o, by = "label") %>%
merge(clay, by = "label") %>%
merge(cec, by = "label")
rownames(hz.data) <- NULL
return(hz.data)
})
# combine horizon data together
hz.data.all <- do.call('rbind', res)
# calculate top and bottom depths from label
spc <- separate(hz.data.all, label, sep = "-", into = c("top", "bottom"))
###
### spc is a tibble with all your data in it
###a
# from here, you can do your analysis e.g. with aqp
# install if needed:
# install.packages(aqp)
# remotes::install_github("ncss-tech/aqp", dependencies=FALSE)
library(aqp)
# promote horizon data to SoilProfileCollection
depths(spc) <- id ~ top + bottom
# plot median/50th percentile SOC for spc truncated to [0,50]
plot(trunc(spc, 0, 50), color = "socQ50")
# plot low/5th percentile
plot(trunc(spc, 0, 50), color = "socQ05")
# plot high/95th percentile
plot(trunc(spc, 0, 50), color = "socQ95")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.