Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Comparision of NASIS or KSSL point data versus SoilGrids point query versus mass-preserving splines using new methods soilDB::fetchSoilGrids and aqp::spc2mpspline
# get latest aqp (1.24) and soilDB (2.5.7) from GitHub
# remotes::install_github("ncss-tech/aqp", dependencies = FALSE)
# remotes::install_github("ncss-tech/soilDB", dependencies = FALSE)
# load packages
library(aqp)
library(soilDB)
# get the classic loafercreek dataset
data(loafercreek)
# or use KSSL
# loafercreek <- fetchKSSL(series = "Loafercreek")
# get the location information
locs <- data.frame(id = profile_id(loafercreek),
site(loafercreek)[,c("y","x")],
stringsAsFactors = FALSE)
# fetch the first 6 pedon locations from SoilGrids
sgr <- fetchSoilGrids(head(locs[complete.cases(locs),]), c("id", "y", "x"))
# get the first 6 pedons from NASIS data (loafercreek dataset)
lcr <- head(filter(loafercreek, complete.cases(locs)))
profile_id(lcr) <- paste0(profile_id(lcr), "_OG")
# promote loafercreek to same spatial CRS
coordinates(lcr) <- ~ x + y
proj4string(lcr) <- "+proj=longlat +datum=WGS84"
# create unique profile ID prior to union
profile_id(sgr) <- paste0(profile_id(sgr),"_SG")
# convert mean clay from SoilGrids to same scale as field clay content
sgr$claycombined <- sgr$clay.mean / 10
lcr$claycombined <- lcr$clay
# combine
lcsg <- aqp::union(list(sgr, lcr))
# compare in horizon-aggregated form
plot(lcsg, color = "claycombined")
# iterating by profile outside of mpspline ensures full depth per profile is captured
lcsg_spline <- aqp::union(profileApply(lcsg, function(x) {
spc2mpspline(x, var_name = "claycombined")
}))
profile_id(lcsg_spline) <- paste0(profile_id(lcsg_spline), "_spline")
# create a combined variable for displaying field, SoilGrids and splines together
lcsg$finalclay <- lcsg$claycombined
lcsg_spline$finalclay <- lcsg_spline$claycombined_spline
# compare in spline'd form
lcsg_all <- aqp::union(list(lcsg, lcsg_spline))
plot(lcsg_all, color = "finalclay", divide.hz = FALSE)
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.