Created
August 20, 2020 21:18
Star
You must be signed in to star a gist
Comparision of NASIS or KSSL point data versus SoilGrids point query versus mass-preserving splines using new methods soilDB::fetchSoilGrids and aqp::spc2mpspline
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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