In ncss-tech/soilReports: R Markdown Reports and Convenience Functions for Soil Survey

## debugging
# params <- list(musym = '3145',
#                cache_file = 'CA792.rda')

# chunk options
knitr::opts_chunk$set(
  message = FALSE,
  warning = FALSE,
  background = '#F7F7F7',
  fig.align = 'center',
  fig.retina = 2,
  dev = 'png',
  tidy = FALSE,
  verbose = FALSE,
  progress = FALSE,
  echo = FALSE
)

library(aqp, quietly = TRUE)
library(soilDB, quietly = TRUE)
library(sharpshootR, quietly = TRUE)
library(latticeExtra, quietly = TRUE)
library(reshape2, quietly = TRUE)
library(tactile, quietly = TRUE)
library(ggplot2, quietly = TRUE)
library(cluster, quietly = TRUE)

# local functions
source('custom.R')

# local configuration
source('config.R')

## re-make cached data
if (is.null(params$cache_file) || !file.exists(params$cache_file)) {
  source('cache-data.R')
} else {
  # load cached data
  load(params$cache_file)
}

## subset pieces

# component + mu records (SPC)
co <- subset(co, musym == params$musym)

# component month (DF)
cm <- subset(cm, subset = coiid %in% profile_id(co))

# component pedon linkage (DF)
cp <- subset(cp, subset = coiid %in% profile_id(co))

# pedons (SPC)
p <- subset(p, profile_id(p) %in% cp$peiid)

# subset extended OSD data here
osds.ac <- subset(osds$climate.annual, subset = series %in% toupper(unique(co$compname)))

# just the SPC
osds <- subset(osds$SPC, profile_id(osds$SPC) %in% toupper(unique(co$compname)))

cotx <- subset(cotx, subset = coiid %in% profile_id(co))

geom <- subset(geom, subset = coiid %in% profile_id(co))

pm <- subset(pm, subset = coiid %in% profile_id(co))

## re-level component labels
co$.label <- factor(
  co$.label,
  levels = co$.label[order(co$comppct_r, decreasing = TRUE)]
)

## add component pedon / component data to pedons

## TODO: check to make sure that there is only a single case of each pedon
any(table(cp$peiid) > 1)

# merge subset component pedons into SPC
site(p) <- cp[, c('peiid', 'coiid', 'rvindicator')]

# look-up associated component label
site(p)$.comp_label <- co$.label[match(p$coiid, site(co)$coiid)]

## TODO: sort component names / labels by decreasing component percent
# pedon convenience label
p$.pedon_label <- p$taxonname

# flag RV copedon
idx <- which(p$rvindicator == 1)
p$.pedon_label[idx] <- paste0(p$.pedon_label[idx], '|R')

# flag OSD copedon
idx <- which(p$pedontype == 'OSD pedon')
p$.pedon_label[idx] <- paste0(p$.pedon_label[idx], '|O')


## establish some reasonable figure widths in inches

# copedon figure
copedon.profile.fig.width <- 2 + (length(p) * 1.5)

# component RV thematic sketches
comp.profile.fig.width <- 2 + (length(co) * 1.5)

# OSD dendrogram for all components that are named series
# this figure is only possible if there are more than 1 series
do.osd.dend <- TRUE
osd.fig.width <- 5 + ((length(co) + length(osds)) * 0.7)

# component comparison by profile_compare
do.comp.comparison <- length(co) > 1

r sprintf("%s: %s", co$musym[1], co$muname[1])
report version r .report.version
r format(Sys.time(), "%Y-%m-%d %H:%M")

This report requires loading several related objects into your NASIS Selected Set, including Area, Legend Mapunit, Correlation, Component Pedon, Pedon and Site Observation.

A useful NASIS query that gets all of the necessary objects is NSSC Pangaea: Area/Legend/Mapunit/DMU/Pedon/Site by areasymbol.

This query has detailed instructions for loading necessary data, and includes only representative data map units. Several other variants of this same query can be used to obtain data based on component or pedon information rather than area symbol.

Components

co.summary <- site(co)[order(co$comppct_r, decreasing = TRUE), c('musym', '.label', 'compkind', 'taxclname', 'hydgrp', 'ecosite_name')]

kableExtra::kable_styling(
  knitr::kable(co.summary, row.names = FALSE, format = 'html'),
   full_width = TRUE, font_size = 11
)

Component Parent Material | Landform

pm.summary <- site(co)[, c('.label', 'landform_string', 'pmkind', 'pmorigin')]

kableExtra::kable_styling(
  knitr::kable(pm.summary, row.names = FALSE, format = 'html'),
   full_width = FALSE
)

Component Text Notes

txt <- cotx[which(cotx$textcat == 'GENSOIL'), ]

txt <- merge(
  site(co)[, c('coiid', '.label', 'comppct_r')],
  txt,
  by = 'coiid',
  all.x = TRUE,
  sort = FALSE
)

txt <- txt[order(txt$comppct_r, decreasing = TRUE), ]

kableExtra::kable_styling(
  knitr::kable(txt[, c('.label', 'textentry')], row.names = FALSE, format = 'html'), full_width = FALSE, font_size = 11
)

txt <- cotx[which(cotx$textcat == 'rep pedon'), ]

txt <- merge(
  site(co)[, c('coiid', '.label', 'comppct_r')],
  txt,
  by = 'coiid',
  all.x = TRUE,
  sort = FALSE
)

txt <- txt[order(txt$comppct_r, decreasing = TRUE), ]

# # attempt to read-in as white-space delimited
# # hmm not going to work, without some clever tricks
# nm <- read.table(textConnection(object = txt$textentry[1]), nrows = 1)
# 
# txt.table <- read.table(textConnection(object = txt$textentry[1]), skip = 3)
# names(txt.table) <- nm

for(i in 1:nrow(txt)) {
  cat('<div id="mu_textnote"><pre>')
  cat(txt$.label[i])
  cat('<hr>')
  cat(txt$textentry[i], sep = '\n')
  cat('</prev></div>')
}

try(print(vizAnnualClimate(osds.ac)$fig))

Component Climate

# comp climate data
nm <- siteNames(co)[grep('^ffd|^maat|^map|^slope|^elev', siteNames(co))]
co.climate <- site(co)[, c('.label', nm)]

# melt
co.climate.long_r <- melt(co.climate, id.vars = '.label', measure.vars = nm[grep('_r', nm)])
co.climate.long_l <- melt(co.climate, id.vars = '.label', measure.vars = nm[grep('_l', nm)])
co.climate.long_h <- melt(co.climate, id.vars = '.label', measure.vars = nm[grep('_h', nm)])

# fix names
names(co.climate.long_r)[3] <- 'rv'

# remove _r
co.climate.long_r$variable <- gsub(pattern = '_r', replacement = '', x = co.climate.long_r$variable, fixed = TRUE)

# row-order is preserved, copy over low / high
co.climate.long_r$low <- co.climate.long_l$value
co.climate.long_r$high <- co.climate.long_h$value

# factor levels / units
co.climate.long_r$variable <- factor(
  co.climate.long_r$variable, 
  levels = c('elev', 'ffd', 'maat', 'map', 'slope'),
  labels = c('Elevation (m)', 'Frost-Free Days', 'MAAT (deg C)', 'MAP (mm)', 'Slope (%)'))

# combine components with area-wide ranges
if(.include.local.ranges) {
  .combined <- rbind(co.climate.long_r, .local_range)
  .combined$.label <- factor(
    .combined$.label, 
    levels = c(levels(co.climate.long_l$.label), unique(.local_range$.label))
  )
}


segplot(
  factor(.label) ~ low + high | variable, 
  data = .combined, 
  centers = rv, 
  draw.bands = FALSE, 
  scales = list(x = list(relation = 'free', rot = 45)),
  strip = strip.custom(bg = grey(0.85)),
  lwd = 2, pch = 15, cex = 1.25,
  par.settings = tactile.theme(),
  layout = c(3, 2),
  as.table = TRUE,
  panel = function(...) {
    panel.grid(h = -1, v = -1)
    panel.segplot(...)
  }
  )

Component | OSD Taxonomic Comparison

# combine OSDs + comp

co.tmp <- co
site(co.tmp)$soilorder <- co.tmp$taxorder
site(co.tmp)$suborder <- co.tmp$taxsuborder
site(co.tmp)$greatgroup <- co.tmp$taxgrtgroup
site(co.tmp)$subgroup <- co.tmp$taxsubgrp

profile_id(co.tmp) <- as.character(co.tmp$.label)
co.tmp <- combine(co.tmp, osds)

# wrapping in try(): there are some cases where there aren't enough data
#                    to compute distance matrix
try(
  SoilTaxonomyDendrogram(co.tmp, cex.taxon.labels = 0.8, width = 0.25, name.style = 'center-center', hz.depths = TRUE, plot.depth.axis = FALSE), 
  silent = TRUE
)

# SoilTaxonomyDendrogram(osds, cex.taxon.labels = 0.8, width = 0.25, name.style = 'center-center', hz.depths = TRUE, plot.depth.axis = FALSE)

Component Diagnostic Features

d.wide <- soilDB:::.diagHzLongtoWide(diagnostic_hz(co), feature = 'featkind', id = 'coiid')
co.tmp <- co
site(co.tmp) <- d.wide

## TODO: something wrong when all data share the same features, likely the re-ordering of vars

try(
  diagnosticPropertyPlot2(co.tmp, v = names(d.wide)[-1], k = 3, grid.label = '.label')
)

Component Comparison

try({
  suppressMessages(
    d <- NCSP(
      co, 
      vars = c('sandtotal_r', 'claytotal_r', 'fragvoltot_r', 'ph1to1h2o_r', 'om_r'),
      # max_d = max(co, v = 'claytotal_r'),
      k = 0,
      # rescale.result = TRUE
    )
  )

  h <- as.hclust(diana(d))

  par(mar = c(0, 0, 0, 0))
  plotProfileDendrogram(co, clust = h, width = 0.25, name.style = 'center-center', hz.depths = TRUE, plot.depth.axis = FALSE, label = '.label', y.offset = 0.15, scaling.factor = 0.007, color = 'claytotal_r')
})

textureTriangleSummary(data.frame(SAND=co$sandtotal_r, SILT=co$silttotal_r, CLAY = co$claytotal_r))

Component Pedons

try({
  par(mar = c(0, 0, 1, 0))
  # component convenience labels
  groupedProfilePlot(
    p,
    groups = '.comp_label',
    label = '.pedon_label',
    group.name.offset = c(-15, -5),
    id.style = 'side',
    group.name.cex = 0.75,
    name.style = 'center-center',
    hz.depths = TRUE,
    plot.depth.axis = FALSE,
    cex.names = 0.66
  )
})

Component RV Evaluation

par(mar=c(0.25, 0.5, 4, 0))

plotSPC(co, label='.label', group.name.cex = 0.75, color='texture', col.label='Texture', col.legend.cex=0.75, name.style = 'center-center', hz.depths = TRUE, plot.depth.axis = FALSE, cex.names = 0.66, width = 0.25)

plotSPC(co, label='.label', group.name.cex = 0.75, color='ksat_r', col.label='Ksat (\U00B5m/sec)', col.legend.cex=0.75, name.style = 'center-center', hz.depths = TRUE, plot.depth.axis = FALSE, cex.names = 0.66, width = 0.25)

plotSPC(co, label='.label', group.name.cex = 0.75, color='om_r', col.label='Organic Matter (%)', col.legend.cex=0.75, name.style = 'center-center', hz.depths = TRUE, plot.depth.axis = FALSE, cex.names = 0.66, width = 0.25)

plotSPC(co, label='.label', group.name.cex = 0.75, color='lep_r', col.label='LEP', col.legend.cex=0.75, name.style = 'center-center', hz.depths = TRUE, plot.depth.axis = FALSE, cex.names = 0.66, width = 0.25)

plotSPC(co, label='.label', group.name.cex = 0.75, color='cec7_r', col.label='CEC @ pH 7 (cmol[+] / kg)', col.legend.cex=0.75, name.style = 'center-center', hz.depths = TRUE, plot.depth.axis = FALSE, cex.names = 0.66, width = 0.25)

plotSPC(co, label='.label', group.name.cex = 0.75, color='ecec_r', col.label='ECEC (cmol[+] / kg)', col.legend.cex=0.75, name.style = 'center-center', hz.depths = TRUE, plot.depth.axis = FALSE, cex.names = 0.66, width = 0.25)

plotSPC(co, label='.label', group.name.cex = 0.75, color='sumbases_r', col.label='Sum of Bases (cmol[+] / kg)', col.legend.cex=0.75, name.style = 'center-center', hz.depths = TRUE, plot.depth.axis = FALSE, cex.names = 0.66, width = 0.25)

Component Low/RV/High Evaluation

# iterate over component IDs
# decreasing component percentage order
co.ids <- site(co)$coiid[order(co$comppct_r, decreasing = TRUE)]
for(comp in co.ids){

  ## subset components: component record ID
  co.sub <- subset(co, site(co)$coiid == comp)

  ## subset pedons: component record ID
  p.sub <- subset(p, coiid == comp)

  ## OSDs: component name
  co.name <- co$compname[match(comp, site(co)$coiid)]
  osds.sub <- subset(osds, id == toupper(co.name))

  # there may not be an OSD yet
  if(length(osds.sub) == 0) {
    # use filler based on deepest component subset
    osds.sub <- emptySPC(co.sub[1, ], top = 0, bottom = max(co.sub)) 

  }

  # there may be no component pedons
  if(length(p.sub) == 0) {
    # use filler based on deepest component subset
    p.sub <- emptySPC(co.sub[1, ], top = 0, bottom = max(co.sub)) 
  }


## overview sketches
OverviewSketches(osds.sub, co.sub, p.sub)

## thematic sketches

# v.co = component basename
# v.p = pedon name

thematicSketches(v.co = 'claytotal', v.p = 'clay', fig.title = 'Clay Content (%)', osds.sub, co.sub, p.sub)

thematicSketches(v.co = 'sandtotal', v.p = 'sand', fig.title = 'Sand Content (%)', osds.sub, co.sub, p.sub)

thematicSketches(v.co = 'fragvoltot', v.p = 'total_frags_pct', fig.title = 'Total Fragments (%)', osds.sub, co.sub, p.sub)

thematicSketches(v.co = 'ph1to1h2o', v.p = 'phfield', fig.title = 'pH 1:1 H2O', osds.sub, co.sub, p.sub)


} # end iterating over components

Component Month

ggplot(cm, aes(month, .label, flodfreqcl)) + geom_tile(aes(fill = flodfreqcl), color='white', lwd=1) + scale_fill_brewer(palette = "Spectral", drop=FALSE, na.value='grey80', name='Flooding Frequency') + scale_x_discrete(drop=FALSE) + theme_minimal() + xlab('') + ylab('') + theme(legend.position="bottom") + guides(fill=guide_legend(nrow=2,byrow=TRUE))

ggplot(cm, aes(month, .label, floddurcl)) + geom_tile(aes(fill = floddurcl), color='white', lwd=1) + scale_fill_brewer(palette = "Spectral", drop=FALSE, na.value='grey80', name='Flooding Duration') + scale_x_discrete(drop=FALSE) + theme_minimal() + xlab('') + ylab('') + theme(legend.position="bottom") + guides(fill=guide_legend(nrow=2,byrow=TRUE))

ggplot(cm, aes(month, .label, pondfreqcl)) + geom_tile(aes(fill = pondfreqcl), color='white', lwd=1) + scale_fill_brewer(palette = "Spectral", drop=FALSE, na.value='grey80', name='Ponding  Frequency') + scale_x_discrete(drop=FALSE) + theme_minimal() + xlab('') + ylab('') + theme(legend.position="bottom") + guides(fill=guide_legend(nrow=2,byrow=TRUE))

ggplot(cm, aes(month, .label, ponddurcl)) + geom_tile(aes(fill = ponddurcl), color='white', lwd=1) + scale_fill_brewer(palette = "Spectral", drop=FALSE, na.value='grey80', name='Ponding  Duration') + scale_x_discrete(drop=FALSE) + theme_minimal() + xlab('') + ylab('') + theme(legend.position="bottom") + guides(fill=guide_legend(nrow=2,byrow=TRUE))