if (FALSE) {
library(aqp)
your.points <- data.frame(id = c("A", "B"),
lat = c(37.9, 38.1),
lon = c(-120.3, -121.5),
stringsAsFactors = FALSE)
x <- try(fetchSoilGrids(your.points))
if (!inherits(x, 'try-error'))
aqp::plotSPC(x, name = NA, color = "socQ50")
# organic carbon stocks use 0-30cm interval
y <- try(fetchSoilGrids(your.points[1, ],
depth_interval = c("0-5", "0-30", "5-15", "15-30"),
variables = c("soc", "bdod", "ocd", "ocs")))
# extract horizons from a SoilProfileCollection where horizon 2 overlaps 1, 3, and 4
h <- aqp::horizons(y)
# "ocs" (organic carbon stock 0-30cm interval)
h[2, ]
h$thickness_meters <- ((h$hzdepb - h$hzdept) / 100)
# estimate "ocs" from modeled organic carbon and bulk density in 0-5, 5-15, 15-30 intervals
# (sum the product of soc, bdod, and thickness in meters)
# (1 gram per cubic decimeter = 1 kilogram per cubic meter)
sum(h$socmean * h$bdodmean * h$thickness_meters, na.rm = TRUE)
# estimate "ocs" from modeled organic carbon density in 0-5, 5-15, 15-30 intervals
# (sum the product of "ocd" and thickness in meters)
sum(h$ocdmean * h$thickness_meters, na.rm = TRUE)
}
Run the code above in your browser using DataLab