Load packages

library(aqp)
library(soilDB)
library(sharpshootR)
library(dendextend)
library(SoilTaxonomy)
library(lattice)
library(kableExtra)

Visualizing soil horizons for each soil series

Based on Official Series Description (OSD)

#Data extracted from Soil Series Descriptions (OSDs)
#Creates a SoilProfileCollection object
soils <- c('Saco', 'Merrimac','Ridgebury','Woodbridge', 'Paxton', 'Hinckley')
osd <- fetchOSD(soils) 


#Dataframe by horizon
#horizons<- as.data.frame(osd@horizons)

# Three different styles
par(mar = c(0,0,1,1), xpd = NA)
plotSPC(osd[1:6, ], cex.names = 1)

par(mar = c(0,0,1,1))
plotSPC(osd, cex.names = 1, name.style = 'center-center', width = 0.3, plot.depth.axis = FALSE, hz.depths = TRUE, fixLabelCollisions = TRUE, hz.depths.offset = 0.06)

par(mar = c(0,0,0,0))
plotSPC(osd, name = 'hzname', plot.order = order(osd$soilseriesiid), 
        cex.names = 1, axis.line.offset = -5, width = 0.33, 
        id.style = 'top', name.style = 'center-center', 
        plot.depth.axis = FALSE, hz.depths = TRUE)

Soils described by the proportion of hillslope positions they are found in

#Get a set of extended site information
s <- fetchOSD(soils, extended = TRUE) 

res <- vizHillslopePosition(s$hillpos)
print(res$fig)

Taxonomically, how do series relate?

osd$soil.depth <- profileApply(osd, estimateSoilDepth, name = 'hzname')

d <- SoilTaxonomyDendrogram(osd, cex.taxon.labels = 0.8, width = 0.33, y.offset = 0.4, shrink = TRUE)

Central tendancy & spread of properties with depth

Estimated lab data for sampled pedons

#Lab data from the Kellogg Soil Survey Laboratory (KSSL)
kssl.data <- fetchKSSL(series = soils)

#Fix capitalization issues
for(i in soils) {
  kssl.data$taxonname[grep(i, kssl.data$taxonname, ignore.case = TRUE)] <- toupper(i)
}


#How many pedons do we have data for?
#NOTE: no lab data for Saco
kableExtra::kable_styling(knitr::kable(table(kssl.data$taxonname), format = 'html'), full_width = FALSE)


g.slab <- slab(kssl.data, taxonname ~ clay + sand + oc + db_od + Ks)
levels(g.slab$variable) <- c('Clay (%)', 'Sand (%)','Organic Carbon (%)', 'Bulk Density', 'Ksat (ROSETTA)') #Set headers
tps <- list(superpose.line=list(col=c('RoyalBlue'), lwd=2))


xyplot(top ~ p.q50 | variable * taxonname, data=g.slab, ylab='Depth', sub='source: KSSL',
       xlab='median bounded by 25th and 75th percentiles',
       lower=g.slab$p.q25, upper=g.slab$p.q75, ylim=c(205,-5),
       panel=panel.depth_function, alpha=0.25, sync.colors=TRUE,
       prepanel=prepanel.depth_function,
       cf=g.slab$contributing_fraction,
       par.strip.text=list(cex=0.8),
       strip=strip.custom(bg=grey(0.85)),
        scales=list(x=list(alternating=1, relation='free'), y=list(alternating=3)),
       par.settings=tps,
       auto.key=list(columns=3, lines=TRUE, points=FALSE, as.table = TRUE)
)

Filling in OSD profile with KSSL data

First- lets look at what the Ksat data (calculated with Rosetta) for pedons analyzed in the KSSL lab look like. Multiple pedons for each soil series were analyzed, so we need to find a way to aggregate multiple profiles/horizons to represent a single series illustration.

library(tidyverse)

data <- fetchKSSL(series = 'Paxton')

df <- as.data.frame(data@horizons)

#Remove profiles with no Ksat data
#NOTE: There is a way to fill in missing values, but for the sake of this exercise I am not.
na <- df %>%
  group_by(pedon_key) %>%
  dplyr::summarize(count_na = sum(!is.na(Ks)))%>%
  filter(count_na > 0)

df <- filter(df, df$pedon_key %in% na$pedon_key)

df <- df %>%
  rename("id" = "pedon_key",
         "top" = "hzn_top",
         "bottom" = "hzn_bot") %>%
  filter(id != "8454") #weird data

#Converting to SoilProfileCollection object
depths(df) <- id ~ top + bottom

par(mar=c(0,0,3,0))
plotSPC(df, name='hzn_desgn', color='Ks', col.palette=viridis::viridis(10), col.label='Paxton Ksat: Described pedons with data')

Horizon depths with Ksat measurements

The profile given in the OSD is a generalized model of what we know to be true on the ground. If we want to use the OSD to represent properties with depth we need to figure out if the depths of the horizons actually sampled can transfer. So let’s inspect horizon mid-points

### do depths of horizons represent the horizons in the OSD?

pax <- fetchOSD('Paxton') 
pax <- as.data.frame(pax@horizons)
pax <- pax %>%
  rename("hzn_desgn" = "hzname")

unique(pax$hzn_desgn)

df2 <- as.data.frame(df@horizons)
unique(df2$hzn_desgn) #list of unique horizon names

#Mid-point depth of any horizon
df2$mid <- (df2$top + df2$bottom)/2
pax$mid <- (pax$top + pax$bottom)/2

df2$type = "KSSL"
pax$type = "OSD"

unique(df2$hzn_desgn)

df2 <- df2 %>%
  filter(hzn_desgn %in% c("Ap", "Bw", "Bw1", "Bw2", "Cd", "Cd1", "Cd2", 'Cd3'))


#To get n= at the bottom of the plot
n <- count(df2, hzn_desgn)
 
ggplot(df2, aes(hzn_desgn, mid, fill = "lightgray")) +
  geom_point(data = pax, aes(hzn_desgn, mid,  fill = type), color = "blue", size = 4)+
  geom_boxplot()+
  theme_classic()+
  geom_text(data = n, aes(y = 0, label = n), color = "darkgray")+
  scale_fill_manual(values = c("lightgray", "blue"), labels = c("KSSL data", "OSD description"))+
  labs(fill = " ",
       title= "Similar horizons with comparable mid-depths  to OSD")+
  ylab("horizon mid-point depth (cm)")+
  xlab("Horizon")

UConn Ksat Soil Series

Amanda Pennino

2023-09-27