## this example investigates the concept of a "median profile"

# required packages
if (require(ape) & require(cluster)) {
  data(sp3)

  # generate a RGB version of soil colors
  # and convert to HSV for aggregation
  sp3$h <- NA
  sp3$s <- NA
  sp3$v <- NA
  sp3.rgb <- with(sp3, munsell2rgb(hue, value, chroma, return_triplets = TRUE))

  sp3[, c('h', 's', 'v')] <- t(with(sp3.rgb, rgb2hsv(r, g, b, maxColorValue = 1)))

  # promote to SoilProfileCollection
  depths(sp3) <- id ~ top + bottom

  # aggregate across entire collection
  a <- slab(sp3, fm = ~ clay + cec + ph + h + s + v, slab.structure = 10)

  # check
  str(a)

  # convert back to wide format
  library(data.table)

  a.wide.q25 <- dcast(a, top + bottom ~ variable, value.var = c('p.q25'))
  a.wide.q50 <- dcast(a, top + bottom ~ variable, value.var = c('p.q50'))
  a.wide.q75 <- dcast(a, top + bottom ~ variable, value.var = c('p.q75'))

  # add a new id for the 25th, 50th, and 75th percentile pedons
  a.wide.q25$id <- 'Q25'
  a.wide.q50$id <- 'Q50'
  a.wide.q75$id <- 'Q75'

  # combine original data with "mean profile"
  vars <- c('top', 'bottom', 'id', 'clay', 'cec', 'ph', 'h', 's', 'v')
  # make data.frame version of sp3
  sp3.df <- as(sp3, 'data.frame')
  sp3.grouped <- rbind(sp3.df[, vars], a.wide.q25[, vars], a.wide.q50[, vars], a.wide.q75[, vars])

  # re-constitute the soil color from HSV triplets
  # convert HSV back to standard R colors
  sp3.grouped$soil_color <- with(sp3.grouped, hsv(h, s, v))

  # give each horizon a name
  sp3.grouped$name <- paste(
    round(sp3.grouped$clay),
    '/' ,
    round(sp3.grouped$cec),
    '/',
    round(sp3.grouped$ph, 1)
  )

  ## perform comparison, and convert to phylo class object
  ## D is rescaled to [0,]
  d <- profile_compare(
                        sp3.grouped,
                        vars = c('clay', 'cec', 'ph'),
                        max_d = 100,
                        k = 0.01,
                        replace_na = TRUE,
                        add_soil_flag = TRUE,
                        rescale.result = TRUE
                      )

  h <- agnes(d, method = 'ward')
  p <- ladderize(as.phylo(as.hclust(h)))

  # look at distance plot-- just the median profile
  plot_distance_graph(d, 12)

  # similarity relative to median profile (profile #12)
  round(1 - (as.matrix(d)[12,] / max(as.matrix(d)[12,])), 2)

  ## make dendrogram + soil profiles
  # first promote to SoilProfileCollection
  depths(sp3.grouped) <- id ~ top + bottom

  # setup plot: note that D has a scale of [0,1]
  par(mar = c(1, 1, 1, 1))

  # get the last plot geometry
  lastPP <- get("last_plot.phylo", envir = .PlotPhyloEnv)

  # the original labels, and new (indexed) order of pedons in dendrogram
  d.labels <- attr(d, 'Labels')

  new_order <- sapply(1:lastPP$Ntip,
                      function(i)
                        which(as.integer(lastPP$xx[1:lastPP$Ntip]) == i))

  # plot the profiles, in the ordering defined by the dendrogram
  # with a couple fudge factors to make them fit
  plot(
    sp3.grouped,
    color = "soil_color",
    plot.order = new_order,
    scaling.factor = 0.01,
    width = 0.1,
    cex.names = 0.5,
    y.offset = max(lastPP$yy) + 0.1,
    add = TRUE
  )
}