Open dylanbeaudette opened 4 years ago
Re-visiting the 2021 data, and comparing with SSURGO / KSSL. Look OK to me.
library(aqp)
library(soilDB)
# KSSL
x <- fetchKSSL(series = 'tulare')
# SSURGO
y <- fetchSDA("compname = 'tulare'")
# grouping ID for later
site(x)$origin <- 'KSSL'
site(y)$origin <- 'SSURGO'
# normalize select hz properties
x$db_thirdbar <- x$db_13b
y$db_thirdbar <- y$dbthirdbar_r
x$om <- x$estimated_om
y$om <- y$om_r
x$name <- x$hzn_desgn
y$name <- y$hzname
# combine
z <- aqp::combine(x, y)
# no fragments
horizons(z)$soil_fraction <- 1
# horizon SOM
hd <- horizonDepths(z)
z$SOC_kg_sq_m <- profileApply(z, FUN = function(i) {
# horizon thickness
i$thick <- i[[hd[2]]] - i[[hd[1]]]
# horizon SOM kg / sq. m.
# horizon thickness (cm) * OM fraction (0-1) * soil fraction (0-1) * Db 1/3 bar (g/cm^3) * unit conversion
res <- i$thick * i$om / 100.0 * i$soil_fraction * i$db_thirdbar * 10
return(res)
})
# sum
z$sum_SOC_kg_sq_m <- profileApply(z, function(i) sum(i$SOC_kg_sq_m, na.rm=TRUE))
# simple plot
par(mar=c(0,0,3,1))
groupedProfilePlot(z, groups = 'origin', color = 'om', group.name.offset = -20, print.id=TRUE)
groupedProfilePlot(z, groups = 'origin', color = 'db_thirdbar', group.name.offset = -20, print.id=TRUE)
groupedProfilePlot(z, groups = 'origin', color = 'SOC_kg_sq_m', group.name.offset = -20, print.id=TRUE)
# summarize: note that one of the SSURGO components is bogus: missing lower horizon data
s <- site(z)[, c('origin', 'sum_SOC_kg_sq_m')]
tapply(s$sum_SOC_kg_sq_m, s$origin, median)
1/3 bar Db values seem way too low.
https://casoilresource.lawr.ucdavis.edu/soil-properties/?prop=db&lat=36.0602&lon=-119.7867&z=9
Checking the current SSURGO data, Tulare components seem too low.
library(aqp)
library(soilDB)
library(sf)
library(terra)
library(viridisLite)
## investigate Db at 1/3 bar
x <- fetchSDA("mukey = '461824'", duplicates = FALSE)
horizons(x)$.hzlab <- sprintf("%s\n%s g/cc", x$hzname, x$dbthirdbar_r)
.a <- list(
label = 'compname',
shrink = TRUE,
scaling.factor = 1.2,
name.style = 'center-center',
cex.names = 0.75
)
options(.aqp.plotSPC.args = .a)
par(mar = c(2, 0, 3, 2))
plotSPC(x, color = 'dbthirdbar_r', name = '.hzlab', col.label = 'Bulk Density at 1/3 Bar (g/cc)')
.txt <- sprintf("National MUSYM: %s", unique(x$nationalmusym))
mtext(.txt, side = 1, at = 6)
The 1/3 bar Db values in KSSL are surprisingly low.
Are these problems, or just the unique landscape position and thick lacustrine deposits?