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

## load packages
library(aqp)
library(soilDB)
library(reshape2)
library(plyr)
library(xtable)
library(Hmisc)
library(latticeExtra)
library(gridExtra)
library(RColorBrewer)
library(sharpshootR)
library(MASS)
library(sf)
library(terra)

## local functions
source('custom.R')

## load configuration
source('config.R')

## load genhz patterns
source(GENHZ_RULES)

## report formatting:
knitr::opts_chunk$set(
  message = FALSE,
  warning = FALSE,
  background = '#F7F7F7',
  out.width = "100%",
  fig.align = 'center',
  dpi = 100,
  fig.retina = 2,
  dev = 'png',
  antialias = 'cleartype',
  tidy = FALSE,
  verbose = FALSE
)

options(width = 100,
        stringsAsFactors = FALSE,
        cache = TRUE)

# save as global options
options(
  p.low.rv.high = p.low.rv.high,
  q.type = q.type,
  ml.profile.smoothing = ml.profile.smoothing
)

## load all pedons from the selected set, do not apply horizonation check / removal
f <- fetchNASIS(rmHzErrors = FALSE,
                nullFragsAreZero = TRUE,
                mixColors = FALSE,
                fill = TRUE)

# ## test with loafercreek
# data(loafercreek)
# f <- loafercreek
# f$longstddecimaldegrees <- f$longstddecimaldegrees
# f$latstddecimaldegrees <- f$latstddecimaldegrees

# init coordinates
initSpatial(f, "EPSG:4326") <- ~ longstddecimaldegrees + latstddecimaldegrees

sf::sf_use_s2(FALSE)

## map unit data: load the official version
if (!is.null(SPATIAL_DSN) && 
    !is.null(SPATIAL_LAYER) && 
    !SPATIAL_DSN %in% c('SSURGO','STATSGO')) {
  mu <- sf::st_read(SPATIAL_DSN, SPATIAL_LAYER)
} else {
  if (is.null(SPATIAL_DSN))
    SPATIAL_DSN <- "SSURGO"

  if (is.null(SPATIAL_LAYER))
    SPATIAL_LAYER <- "mupolygon" 

  mukey <- soilDB::SDA_spatialQuery(as(f, 'sf'), 
                                    what = SPATIAL_LAYER, 
                                    db = SPATIAL_DSN)
  nmusym <- soilDB::SDA_query(paste0("SELECT DISTINCT mukey, nationalmusym, musym, areasymbol FROM mapunit 
                                      INNER JOIN legend ON mapunit.lkey = legend.lkey WHERE mukey IN ",
                              soilDB::format_SQL_in_statement(mukey$mukey)))
  mu <- merge(mukey, nmusym, by = "mukey", sort = FALSE)

  ## FULL EXTENT OF nationalmusym
  # mu <- sf::st_as_sf(soilDB::fetchSDA_spatial(
  #   nmusym$nationalmusym,
  #   by.col = "nationalmusym",
  #   verbose = FALSE,
  #   db = SPATIAL_DSN,
  #   add.fields = c("legend.areasymbol", "mapunit.musym")
  # ))

  mu$MUSYM <- mu$musym
  mu$AREASYMBOL <- mu$areasymbol 
}

# transform from GCS to CRS of map unit linework
f_sp <- st_transform(as(f, 'sf'), st_crs(mu))

## overlay with map unit polys, and clean-up
f$musym <- mu$MUSYM[sapply(st_intersects(f_sp, mu), function(x)
  if (length(x) > 0) {
    x
  } else{
    NA
  })]

SUBSET_RULE <- tolower(SUBSET_RULE)
if (SUBSET_RULE == 'pedon.id.list') {
  subset.idx <- which(f$upedonid %in% pedon_set)
}

## generate index to subset using regular expression
if (SUBSET_RULE == 'pattern') {
  subset.idx <- grep(pattern = TAXONNAME_PATTERN, f$taxonname, ignore.case = TRUE)
}

if (SUBSET_RULE == 'musymtaxon') {
  subset.idx <- grepl(pattern = TAXONNAME_PATTERN, f$taxonname, ignore.case = TRUE) &
                grepl(pattern = TAXONKIND_PATTERN, f$taxonkind, ignore.case = TRUE) &
                grepl(pattern = MUSYM_PATTERN, f$musym, ignore.case = TRUE)
}

if (SUBSET_RULE == 'musym') {
  subset.idx <- grep(pattern = MUSYM_PATTERN, f$musym, ignore.case = TRUE)
}

# perform subset
f <- f[subset.idx, ]

if (length(f) == 0) {
  stop('No pedons after applying subset rules from config.R')
}

# apply genhz rules based on exact match of TAXONNAME_PATTERN
# TODO: this should allow application of patterns to subsets of SPC
if (TAXONNAME_PATTERN %in% names(gen.hz.rules)) {
  p <- gen.hz.rules[[TAXONNAME_PATTERN]]$p
  n <- gen.hz.rules[[TAXONNAME_PATTERN]]$n
  f <- generalizeHz(f, new = n, pattern = p)
}

## update diagnostic feature slot
# join upedonid + additional information into diagnostic table: this is kind of wasteful
f.diagnostic <- join(site(f)[, c('upedonid', 'peiid', 'musym', 'taxonname')], 
                     diagnostic_hz(f), by = 'peiid', type = 'left')

# remove records where diag_kind is NA
missing.featkind <- which(is.na(f.diagnostic$featkind))
if(length(missing.featkind) > 0)
    f.diagnostic <- f.diagnostic[-missing.featkind, ]

# copy diagnostic data into @diagnostic as list
diagnostic_hz(f) <- f.diagnostic

## overlay point data with raster data
# extract site+coordinates for overlay
f.sp <- as(f, 'sf')

# iterate over rasters, and extract values at pedon locations
r <- lapply(RASTER_LIST, terra::rast)
l.res <- lapply(lapply(r, function(x) {
  terra::extract(x, terra::vect(sf::st_transform(f.sp, sf::st_crs(x))))
}), `[[`, 2)

# convert to DF
l.res <- as.data.frame(l.res, stringsAsFactors=FALSE)

# order is preserved so we can include peiid from sites
l.res$peiid <- f.sp$peiid

## add sampled GIS data to site-level attributes in SPC
site(f) <- l.res

## convert some columns to factors and set levels
if (!is.null(f$gis_geomorphons)) {
  # set geomorphons levels
  f$gis_geomorphons <- factor(
    f$gis_geomorphons,
    levels = 1:10,
    labels = c(
      'flat',
      'summit',
      'ridge',
      'shoulder',
      'spur',
      'slope',
      'hollow',
      'footslope',
      'valley',
      'depression'
    )
  )
}

### TODO: this will cause problems when the rules file isn't up to date
# set GLH levels from original rules
# f$genhz <- factor(f$genhz, levels=gen.hz.rules[[comp]]$n)

# set GHL levels from depths
f$genhz <- factor(f$genhz, levels = guessGenHzLevels(f, 'genhz')$levels)

# compute depth-class information
sdc <- getSoilDepthClass(f)
site(f) <- sdc

### TODO: un-pack this function
# compute summaries
s <- summarize.component(f)

# determine max number of profiles:
max.comp.profiles <- s$n

Component Report

r format(Sys.time(), "%Y-%m-%d")

r paste(MUSYM_PATTERN, TAXONNAME_PATTERN, TAXONKIND_PATTERN)

ranges are (r p.low.rv.high) percentiles

Taxon Names and Pedon Types

Check to make sure that pedons used within this report have been correctly assigned to this component. If not, please fix in NASIS.

wzxhzdk:1

MUSYM Summary

wzxhzdk:2

Hillslope Position Summary

wzxhzdk:3

Geomorphic Component Summaries

wzxhzdk:4

wzxhzdk:5

Geomorphons Summary

this.data <- categorical.prop.table(f$gis_geomorphons)
this.align <- rep('c', times = ncol(this.data) + 1)

print(
  xtable(this.data, align = this.align),
  type = 'html',
  include.rownames = FALSE,
  table.placement = "H",
  caption.placement = "top",
  html.table.attributes = 'cellpadding="3" cellspacing="3"'
)

Drainage Class Summary

wzxhzdk:7

Surface Shape Summary

wzxhzdk:8

Ecosite Summary

wzxhzdk:9

Generalized Horizon Classification

These tables describe the mapping between field-described horizonation (top row) and generalized horizonation (first column). Numbers describe the number of times any given field-described horizon has been allocated to a generalized horizon. If present, values in the "NA" row should be further investigated.

wzxhzdk:10

cols <- c(rev(brewer.pal(11, 'Spectral')))
col.palette <- colorRampPalette(cols)

this.data <- t(s$ct)
this.data[this.data == 0] <- NA
if (ncol(this.data) > 0) {
  levelplot(
    this.data,
    col.regions = col.palette,
    colorkey = list(tick.number = 15),
    xlab = 'Original Horizon Designation',
    ylab = 'GHL',
    main = 'GHL Assignment Evaluation',
    scales = list(alternating = 3),
    panel = function(x, y, z, ...) {
      panel.levelplot(x, y, z, ...)
      idx <- which(!is.na(z))
      panel.text(x[idx], y[idx], z[idx], font = 2)
      panel.abline(h = seq(from = 0.5, to = length(y), by = 1),
                   col = grey(0.45))
      panel.abline(v = seq(from = 0.5, to = length(x), by = 1),
                   col = grey(0.45))
    }
  )
}

GHL assignment as a network graph.

this.data <- t(s$ct)
# convert contingency table -> adj. matrix
m <- genhzTableToAdjMat(this.data)
if (any(m > 0)) {
  # plot using a function from the sharpshootR package
  par(mar = c(1, 1, 1, 1))
  plotSoilRelationGraph(
    m,
    graph.mode = 'directed',
    edge.arrow.size = 0.5,
    vertex.label.family = 'sans'
  )
}

# clay box-whisker plot, grouped by genhz, over-printed with original hz names
# subset data
h.i <- horizons(f)
h.i.sub <- subset(h.i, subset = !is.na(clay), drop = TRUE)
# hack: reset factor levels, to accomodate filtered O horizons
h.i.sub$genhz <- factor(h.i.sub$genhz)

# plotting style
tps <- list(
  box.umbrella = list(col = grey(0.4)),
  box.rectangle = list(col = grey(0.4)),
  box.dot = list(col = grey(0.4), cex = 0.75),
  plot.symbol = list(col = grey(0.4), cex = 0.5)
)
# plot
print(bwplot(
  genhz ~ clay,
  data = h.i.sub,
  main = f,
  par.settings = tps
) + layer(
  panel.text(
    x = h.i.sub$clay,
    y = jitter(as.numeric(h.i.sub$genhz), factor = 1.5),
    label = h.i.sub$hzname,
    cex = 0.75,
    font = 2,
    col = 'RoyalBlue'
  )
))

Maximum-Likelihood Horizonation

The figure below describes the most likely horizonation, based on the collection of pedons associated with this component. This is only an estimate, expert knowledge should be used to adjust these values as needed. When pedon numbers are low or horizonation is not consistent, overlap can occur. Values in square brackets are related to Brier Scores, smaller values suggest more consistent horizonation within the collection.

trellis.par.set(list(superpose.line = list(lwd = 2)))
print(s$ml.hz.plot)

Component Profile Plot

These profile sketches represent the entire collection of named components within the selected set, ordered by map unit symbol.

# this resets the default image width according to the number of profiles
knitr::opts_chunk$set(fig.width = max.comp.profiles * 1.25)
knitr::opts_chunk$set(fig.height = 4)

wzxhzdk:16

Texture Class Summary Tables

These tables describe the frequency of textural classes, summarized by component, map unit and generalized horizon. Values within parenthesis are the fraction of horizons associated with each texture class.

wzxhzdk:17

Morphologic Summary Tables

These table describe low-rv-high values for morphologic properties, summarized by component. The low values are the r p.low.rv.high[1] percentile, RV values are the r p.low.rv.high[2] percentile, and the high values are the r p.low.rv.high[3] percentile.

wzxhzdk:18

Aggregate Color Summary, dry

par(mar = c(4.5, 2, 0, 0))
f$genhz <- as.character(f$genhz)
f$genhz[is.na(f$genhz)] <- "<not-used>"
aggregateColorPlot(
  aggregateColor(f, groups = 'genhz', col = 'dry_soil_color'),
  label.font = 2,
  label.cex = 0.95,
  print.n.hz = TRUE
)

Aggregate Color Summary, moist

par(mar = c(4.5, 2, 0, 0))
aggregateColorPlot(
  aggregateColor(f, groups = 'genhz', col = 'moist_soil_color'),
  label.font = 2,
  label.cex = 0.95,
  print.n.hz = TRUE
)

Morphologic Summary by Map Unit

Whiskers extend from the 5th to 95th percentiles, the body represents the 25th through 75th percentiles, and the dot is the 50th percentile.

print(s$pmg)

Surface Fragment Summary Tables

These table describe low-rv-high values for surface rock fragments, summarized by component and map unit. The low values are the r p.low.rv.high[1] percentile, RV values are the r p.low.rv.high[2] percentile, and the high values are the r p.low.rv.high[3] percentile.

this.data <- s$sf
this.align <- c('l', rep('c', times = ncol(this.data)))
print(
  xtable(this.data, align = this.align),
  type = 'html',
  table.placement = "H",
  caption.placement = "top",
  include.rownames = FALSE,
  html.table.attributes = 'cellpadding="3" cellspacing="3"'
)

Diagnostic feature summary

The low values are the r p.low.rv.high[1] percentile, RV values are the r p.low.rv.high[2] percentile, and the high values are the r p.low.rv.high[3] percentile.

this.data <- s$dt
if (!is.null(this.data)) {
  this.align <- c('l', rep('c', times = ncol(this.data)))
  print(
    xtable(this.data, align = this.align),
    type = 'html',
    table.placement = "H",
    caption.placement = "top",
    include.rownames = FALSE,
    html.table.attributes = 'cellpadding="3" cellspacing="5"'
  )
}

diagnosticPropertyPlot2(
  f,
  v = c(
    'lithic.contact',
    'paralithic.contact',
    'argillic.horizon',
    'cambic.horizon',
    'ochric.epipedon',
    'mollic.epipedon',
    'very.shallow',
    'shallow',
    'mod.deep',
    'deep',
    'very.deep'
  ),
  k = 3
)

Pedon GIS Summary

The low values are the r p.low.rv.high[1] percentile, RV values are the r p.low.rv.high[2] percentile, and the high values are the r p.low.rv.high[3] percentile. These values were sampled from raster data sources, at each pedon location. Arrows on the circular histogram of field-measured aspect values are related to percentiles and "mean resultant length", on a circular basis. Grey arrows are the r p.low.rv.high[1] and r p.low.rv.high[3] percentiles and the red arrow is the r p.low.rv.high[2] percentile. Longer arrows suggest an aspect-affected pattern or aspect-biased sampling site selection.

this.data <- s$pg
this.align <- rep('c', times = ncol(this.data) + 1)
i.xt <- xtable(this.data, align = this.align)
digits(i.xt) <- 0
print(
  i.xt,
  type = 'html',
  include.rownames = FALSE,
  table.placement = "H",
  caption.placement = "top",
  html.table.attributes = 'cellpadding="3" cellspacing="3"'
)

# this resets the default image width according to the number of profiles
knitr::opts_chunk$set(fig.width = 4.5)

par(mar = c(0, 0, 0, 0))
aspect.plot(
  f$aspect,
  q = p.low.rv.high,
  plot.title = TAXONNAME_PATTERN,
  pch = 21,
  bg = 'RoyalBlue',
  col = 'black',
  arrow.col = c('grey', 'red', 'grey')
)

try(unlink('this.component.Rda'))

---

This document is based on aqp version r utils::packageDescription("aqp", field="Version") and soilDB version r utils::packageDescription("soilDB", field="Version").

ncss-tech/soilReports documentation built on March 6, 2025, 1:15 a.m.