library(aqp)
library(soilDB)
library(sf)
library(sharpshootR)
library(SoilTaxonomy)


## possible AOIs defined using a bounding-box via SoilWeb

# WI: many ties 
bb <- '-97.0983 39.3808,-97.0983 39.4127,-97.0282 39.4127,-97.0282 39.3808,-97.0983 39.3808'

# KS069
bb <- '-100.5534 37.9177,-100.5534 37.9822,-100.4389 37.9822,-100.4389 37.9177,-100.5534 37.9177'



## assemble AOI polygon into WKT
wkt <- sprintf('POLYGON((%s))', bb)

## init sf polygon
# GCS WGS84
x <- st_as_sfc(wkt, crs = 4326)


## get overlapping map unit keys
# could also use SDA_query() with more elaborate SQL
m <- SDA_spatialQuery(x, what = 'mukey')

## compose SQL to return component details for these map unit keys
# return only:
# * map units overlapping with BBOX
# * major components
# * no misc. areas that might share name with a poorly-named soil series
sql <- sprintf(
  "SELECT mukey, cokey, compname, compkind, comppct_r 
  FROM component 
  WHERE mukey IN %s 
  -- AND majcompflag = 'Yes'
  AND compkind != 'Miscellaneous area'
  ", format_SQL_in_statement(as.integer(m$mukey))
)

## send to SDA, result is a data.frame
s <- SDA_query(sql)



## get OSD morphology + extended summaries 
osd <- fetchOSD(unique(s$compname), extended = TRUE)


## check out results
str(osd, 1)
## List of 18
##  $ SPC             :Formal class 'SoilProfileCollection' [package "aqp"] with 8 slots
##  $ competing       :'data.frame':    87 obs. of  3 variables:
##  $ geog_assoc_soils:'data.frame':    74 obs. of  2 variables:
##  $ geomcomp        :'data.frame':    14 obs. of  9 variables:
##  $ hillpos         :'data.frame':    14 obs. of  8 variables:
##  $ mtnpos          : logi FALSE
##  $ terrace         :'data.frame':    4 obs. of  5 variables:
##  $ flats           :'data.frame':    7 obs. of  7 variables:
##  $ shape_across    :'data.frame':    14 obs. of  8 variables:
##  $ shape_down      :'data.frame':    14 obs. of  8 variables:
##  $ pmkind          :'data.frame':    32 obs. of  5 variables:
##  $ pmorigin        :'data.frame':    5 obs. of  5 variables:
##  $ mlra            :'data.frame':    107 obs. of  4 variables:
##  $ ecoclassid      :'data.frame':    67 obs. of  5 variables:
##  $ climate.annual  :'data.frame':    112 obs. of  12 variables:
##  $ climate.monthly :'data.frame':    336 obs. of  14 variables:
##  $ NCCPI           :'data.frame':    14 obs. of  16 variables:
##  $ soilweb.metadata:'data.frame':    23 obs. of  2 variables:
# single iteration of hydrologic ordering
h1 <- hydOrder(osd$shape_across, g = 'shape', clust = TRUE)

# return trace log for eval of objective function
# increase max iterations
h2 <- iterateHydOrder(osd$shape_across, 'shape', maxIter = 100, verbose = TRUE, trace = TRUE)

# compare
h1$match.rate
## [1] 1
h2$match.rate
## [1] 1
# inspect objective function evolution
tr <- h2$trace

obj <- sapply(tr, '[[', 'obj')
hist(obj, las = 1, xlab = 'Objective Function\nLower is Better')

Histogram of obj Objective Function Lower is Better Frequency 0 20 40 60 80 0 1 2 

match.rate <- sapply(tr, '[[', 'match.rate')
hist(match.rate, las = 1, xlab = 'Matching Rate')

Histogram of match.rate Matching Rate Frequency 0.85 0.90 0.95 1.00 0 1 2

Note that profiles deeper than 180cm are marked as truncated via ragged bottoms.

# provide additional arguments to aqp::plotSPC() via options
options(.aqp.plotSPC.args = list(max.depth = 190))

par(mar = c(1, 0, 1, 2))
plotGeomorphCrossSection(osd, type = 'line', maxIter = 100, j.amount = 0.05, verbose = TRUE)

Ap Bt1 Bt2 Bk1 Bk2 C SPEARVILLE Ap AB Bt1 Bt2 BCk Ck C HARNEY Ap Bt BCk1 BCk2 C RICHFIELD A Bw1 Bw2 Bk DUROC Ap1 Ap2 Bw Bk1 Bk2 Bk3 Bk4 ULYSSES A Bkk Btkk 2Btk1 2Btk2 MANSKER A BA Bw BC C ULY A AC C1 C2 C3 C4 COLY A1 A2 Bk C PENDEN Ap Bw Bk1 Bk2 Bk3 2Bk4 2Bk5 BUFFALO PARK A Bss1 Bss2 BC 2C1 2C2 NESS Ap Bt1 Bt2 Bt3 Bt4 BC C LODGEPOLE Ap BA Bw Bwss1 Bwss2 Bwss3 Bssk1 Bssk2 Bk BCk FETERITA A BA Bt1 Bt2 Bt3 2C PLEASANT 200 cm 175 cm 150 cm 125 cm 100 cm 75 cm 50 cm 25 cm 0 cm 0.0 0.2 0.4 0.6 0.8 1.0 Summit Shoulder Backslope Footslope Toeslope Proportion 

plotGeomorphCrossSection(osd, type = 'bar', maxIter = 100, j.amount = 0.05, verbose = TRUE)

Ap Bt1 Bt2 Bk1 Bk2 C SPEARVILLE Ap AB Bt1 Bt2 BCk Ck C HARNEY Ap Bt BCk1 BCk2 C RICHFIELD A Bw1 Bw2 Bk DUROC Ap1 Ap2 Bw Bk1 Bk2 Bk3 Bk4 ULYSSES A Bkk Btkk 2Btk1 2Btk2 MANSKER A BA Bw BC C ULY A AC C1 C2 C3 C4 COLY A1 A2 Bk C PENDEN Ap Bw Bk1 Bk2 Bk3 2Bk4 2Bk5 BUFFALO PARK A Bss1 Bss2 BC 2C1 2C2 NESS Ap Bt1 Bt2 Bt3 Bt4 BC C LODGEPOLE Ap BA Bw Bwss1 Bwss2 Bwss3 Bssk1 Bssk2 Bk BCk FETERITA A BA Bt1 Bt2 Bt3 2C PLEASANT 200 cm 175 cm 150 cm 125 cm 100 cm 75 cm 50 cm 25 cm 0 cm 12 5 11 3 14 7 13 2 9 1 8 6 4 10 typic argiustolls typic argiustolls aridic argiustolls pachic haplustolls torriorthentic haplustolls calcidic paleustolls typic haplustolls typic ustorthents typic calciustolls aridic haplustepts ustic epiaquerts vertic argiaquolls aridic epiaquerts torrertic argiustolls Summit Shoulder Backslope Footslope Toeslope Proportion 

par(mar = c(1, 0, 1, 2))
plotGeomorphCrossSection(osd, type = 'line', clust = FALSE)

Ap Bt1 Bt2 Bk1 Bk2 C SPEARVILLE Ap AB Bt1 Bt2 BCk Ck C HARNEY Ap Bt BCk1 BCk2 C RICHFIELD Ap1 Ap2 Bw Bk1 Bk2 Bk3 Bk4 ULYSSES A Bkk Btkk 2Btk1 2Btk2 MANSKER A BA Bw BC C ULY A AC C1 C2 C3 C4 COLY A1 A2 Bk C PENDEN Ap Bw Bk1 Bk2 Bk3 2Bk4 2Bk5 BUFFALO PARK A Bw1 Bw2 Bk DUROC A Bss1 Bss2 BC 2C1 2C2 NESS Ap Bt1 Bt2 Bt3 Bt4 BC C LODGEPOLE Ap BA Bw Bwss1 Bwss2 Bwss3 Bssk1 Bssk2 Bk BCk FETERITA A BA Bt1 Bt2 Bt3 2C PLEASANT 175 cm 150 cm 125 cm 100 cm 75 cm 50 cm 25 cm 0 cm 0.0 0.2 0.4 0.6 0.8 1.0 Summit Shoulder Backslope Footslope Toeslope Proportion 

plotGeomorphCrossSection(osd, type = 'bar', clust = FALSE)

Ap Bt1 Bt2 Bk1 Bk2 C SPEARVILLE Ap AB Bt1 Bt2 BCk Ck C HARNEY Ap Bt BCk1 BCk2 C RICHFIELD Ap1 Ap2 Bw Bk1 Bk2 Bk3 Bk4 ULYSSES A Bkk Btkk 2Btk1 2Btk2 MANSKER A BA Bw BC C ULY A AC C1 C2 C3 C4 COLY A1 A2 Bk C PENDEN Ap Bw Bk1 Bk2 Bk3 2Bk4 2Bk5 BUFFALO PARK A Bw1 Bw2 Bk DUROC A Bss1 Bss2 BC 2C1 2C2 NESS Ap Bt1 Bt2 Bt3 Bt4 BC C LODGEPOLE Ap BA Bw Bwss1 Bwss2 Bwss3 Bssk1 Bssk2 Bk BCk FETERITA A BA Bt1 Bt2 Bt3 2C PLEASANT 175 cm 150 cm 125 cm 100 cm 75 cm 50 cm 25 cm 0 cm 12 5 11 14 7 13 2 9 1 3 8 6 4 10 typic argiustolls typic argiustolls aridic argiustolls torriorthentic haplustolls calcidic paleustolls typic haplustolls typic ustorthents typic calciustolls aridic haplustepts pachic haplustolls ustic epiaquerts vertic argiaquolls aridic epiaquerts torrertic argiustolls Summit Shoulder Backslope Footslope Toeslope Proportion 

o <- reconcileOSDGeomorph(osd, 'hillpos')
res <- vizHillslopePosition(o$geom, verbose = TRUE)
print(res$fig)

Proportion SPEARVILLE HARNEY RICHFIELD DUROC ULYSSES MANSKER ULY COLY PENDEN BUFFALO PARK NESS LODGEPOLE FETERITA PLEASANT 0.0 0.2 0.4 0.6 0.8 1.0 13 372 156 172 415 166 657 475 244 74 87 302 84 100 0 1.04 1.51 1.03 1.4 0.78 0.7 0.35 0.04 0 0.48 0.57 0 0 N H Toeslope Footslope Backslope Shoulder Summit 

par(mar = c(1, 0, 1, 1))
idx <- match(hydOrder(o$geom, g = 'hillpos', clust = FALSE), profile_id(o$SPC))
plotSPC(o$SPC, plot.order = idx, width = 0.35, name.style = 'center-center', cex.names = 0.65)

Ap Bt1 Bt2 Bk1 Bk2 C SPEARVILLE Ap AB Bt1 Bt2 BCk Ck C HARNEY Ap Bt BCk1 BCk2 C RICHFIELD Ap1 Ap2 Bw Bk1 Bk2 Bk3 Bk4 ULYSSES A Bkk Btkk 2Btk1 2Btk2 MANSKER A BA Bw BC C ULY A AC C1 C2 C3 C4 COLY A1 A2 Bk C PENDEN Ap Bw Bk1 Bk2 Bk3 2Bk4 2Bk5 BUFFALO PARK A Bw1 Bw2 Bk DUROC A Bss1 Bss2 BC 2C1 2C2 NESS Ap Bt1 Bt2 Bt3 Bt4 BC C LODGEPOLE Ap BA Bw Bwss1 Bwss2 Bwss3 Bssk1 Bssk2 Bk BCk FETERITA A BA Bt1 Bt2 Bt3 2C PLEASANT 175 cm 150 cm 125 cm 100 cm 75 cm 50 cm 25 cm 0 cm 

o <- reconcileOSDGeomorph(osd, 'geomcomp')
res <- vizGeomorphicComponent(o$geom, verbose = TRUE)
print(res$fig)

Proportion SPEARVILLE RICHFIELD ULYSSES BUFFALO PARK MANSKER HARNEY ULY COLY PENDEN DUROC PLEASANT LODGEPOLE FETERITA NESS 0.0 0.2 0.4 0.6 0.8 1.0 13 66 285 74 27 372 638 461 244 134 102 187 82 78 0 1.12 2.09 0.99 1.31 0.93 0.47 0.24 0.04 0.46 0.14 0.17 0 0 N H Interfluve Crest Head Slope Nose Slope Side Slope Base Slope 

par(mar = c(1, 0, 1, 1))
idx <- match(hydOrder(o$geom, g = 'geomcomp', clust = FALSE), profile_id(o$SPC))
plotSPC(o$SPC, plot.order = idx, width = 0.35, name.style = 'center-center', cex.names = 0.65)

Ap Bt1 Bt2 Bk1 Bk2 C SPEARVILLE Ap Bt BCk1 BCk2 C RICHFIELD Ap AB Bt1 Bt2 BCk Ck C HARNEY Ap1 Ap2 Bw Bk1 Bk2 Bk3 Bk4 ULYSSES A Bkk Btkk 2Btk1 2Btk2 MANSKER Ap Bw Bk1 Bk2 Bk3 2Bk4 2Bk5 BUFFALO PARK A BA Bw BC C ULY A AC C1 C2 C3 C4 COLY A1 A2 Bk C PENDEN A Bw1 Bw2 Bk DUROC A BA Bt1 Bt2 Bt3 2C PLEASANT Ap Bt1 Bt2 Bt3 Bt4 BC C LODGEPOLE Ap BA Bw Bwss1 Bwss2 Bwss3 Bssk1 Bssk2 Bk BCk FETERITA A Bss1 Bss2 BC 2C1 2C2 NESS 175 cm 150 cm 125 cm 100 cm 75 cm 50 cm 25 cm 0 cm 

o <- reconcileOSDGeomorph(osd, 'flats')
res <- vizFlatsPosition(o$geom, verbose = TRUE)
print(res$fig)

Proportion MANSKER ULYSSES RICHFIELD FETERITA DUROC LODGEPOLE NESS 0.0 0.2 0.4 0.6 0.8 1.0 89 9 7 4 6 180 21 0.79 0.99 0 0.81 0 0 0 N H Dip Talf Flat Rise 

par(mar = c(1, 0, 1, 2))
idx <- match(hydOrder(o$geom, g = 'flats', clust = FALSE), profile_id(o$SPC))
plotSPC(o$SPC, plot.order = idx, width = 0.35, name.style = 'center-center', cex.names = 0.65)

A Bkk Btkk 2Btk1 2Btk2 MANSKER Ap1 Ap2 Bw Bk1 Bk2 Bk3 Bk4 ULYSSES Ap Bt BCk1 BCk2 C RICHFIELD Ap BA Bw Bwss1 Bwss2 Bwss3 Bssk1 Bssk2 Bk BCk FETERITA A Bw1 Bw2 Bk DUROC Ap Bt1 Bt2 Bt3 Bt4 BC C LODGEPOLE A Bss1 Bss2 BC 2C1 2C2 NESS 175 cm 150 cm 125 cm 100 cm 75 cm 50 cm 25 cm 0 cm 

o <- reconcileOSDGeomorph(osd, 'terrace')
res <- vizTerracePosition(o$geom, verbose = TRUE)
print(res$fig)

Proportion DUROC NESS RICHFIELD ULYSSES 0.0 0.2 0.4 0.6 0.8 1.0 15 4 10 8 0 0 0 0 N H Tread Riser 

par(mar = c(1, 0, 1, 1))
idx <- match(hydOrder(o$geom, g = 'terrace', clust = FALSE), profile_id(o$SPC))
plotSPC(o$SPC, plot.order = idx, width = 0.35, name.style = 'center-center', cex.names = 0.65)

A Bw1 Bw2 Bk DUROC A Bss1 Bss2 BC 2C1 2C2 NESS Ap Bt BCk1 BCk2 C RICHFIELD Ap1 Ap2 Bw Bk1 Bk2 Bk3 Bk4 ULYSSES 175 cm 150 cm 125 cm 100 cm 75 cm 50 cm 25 cm 0 cm 

o <- reconcileOSDGeomorph(osd, 'mtnpos')
res <- vizMountainPosition(o$geom, verbose = TRUE)
print(res$fig)
par(mar = c(1, 0, 1, 0))
idx <- match(hydOrder(o$geom, g = 'mtnpos', clust = FALSE), profile_id(o$SPC))
plotSPC(o$SPC, plot.order = idx, width = 0.35, name.style = 'center-center', cex.names = 0.65)
o <- reconcileOSDGeomorph(osd, 'shape_across')
res <- vizSurfaceShape(o$geom, title = 'Shape Across', verbose = TRUE)
print(res$fig)

Shape Across Proportion HARNEY PENDEN ULY COLY BUFFALO PARK RICHFIELD ULYSSES MANSKER SPEARVILLE DUROC LODGEPOLE NESS PLEASANT FETERITA 0.0 0.2 0.4 0.6 0.8 1.0 373 245 671 760 74 175 426 179 13 212 422 192 192 87 1 0.95 0.82 0.78 0.75 0.89 0.94 0.31 0 0.96 0.43 0.15 0.15 0 N H Concave Linear Convex Complex Undulating 

par(mar = c(1, 0, 1, 1))
idx <- match(hydOrder(o$geom, g = 'shape', clust = FALSE), profile_id(o$SPC))
plotSPC(o$SPC, plot.order = idx, width = 0.35, name.style = 'center-center', cex.names = 0.65)

Ap AB Bt1 Bt2 BCk Ck C HARNEY A1 A2 Bk C PENDEN A BA Bw BC C ULY A AC C1 C2 C3 C4 COLY Ap Bw Bk1 Bk2 Bk3 2Bk4 2Bk5 BUFFALO PARK Ap Bt BCk1 BCk2 C RICHFIELD Ap1 Ap2 Bw Bk1 Bk2 Bk3 Bk4 ULYSSES A Bkk Btkk 2Btk1 2Btk2 MANSKER Ap Bt1 Bt2 Bk1 Bk2 C SPEARVILLE A Bw1 Bw2 Bk DUROC Ap Bt1 Bt2 Bt3 Bt4 BC C LODGEPOLE A Bss1 Bss2 BC 2C1 2C2 NESS A BA Bt1 Bt2 Bt3 2C PLEASANT Ap BA Bw Bwss1 Bwss2 Bwss3 Bssk1 Bssk2 Bk BCk FETERITA 175 cm 150 cm 125 cm 100 cm 75 cm 50 cm 25 cm 0 cm 

o <- reconcileOSDGeomorph(osd, 'shape_down')
res <- vizSurfaceShape(o$geom, title = 'Shape Down', verbose = TRUE)
print(res$fig)

Shape Down Proportion BUFFALO PARK COLY PENDEN HARNEY RICHFIELD ULYSSES ULY SPEARVILLE DUROC MANSKER PLEASANT NESS FETERITA LODGEPOLE 0.0 0.2 0.4 0.6 0.8 1.0 74 760 245 373 175 426 671 13 212 179 192 192 87 422 0.84 1.12 0.99 1 0.9 1.36 1.33 0.62 0.85 1.13 1.09 0.15 0.09 0 N H Concave Linear Convex Complex Undulating 

par(mar = c(1, 0, 1, 1))
idx <- match(hydOrder(o$geom, g = 'shape', clust = FALSE), profile_id(o$SPC))
plotSPC(o$SPC, plot.order = idx, width = 0.35, name.style = 'center-center', cex.names = 0.65)

Ap Bw Bk1 Bk2 Bk3 2Bk4 2Bk5 BUFFALO PARK A AC C1 C2 C3 C4 COLY A1 A2 Bk C PENDEN Ap AB Bt1 Bt2 BCk Ck C HARNEY Ap Bt BCk1 BCk2 C RICHFIELD Ap1 Ap2 Bw Bk1 Bk2 Bk3 Bk4 ULYSSES A BA Bw BC C ULY A Bkk Btkk 2Btk1 2Btk2 MANSKER Ap Bt1 Bt2 Bk1 Bk2 C SPEARVILLE A Bw1 Bw2 Bk DUROC A BA Bt1 Bt2 Bt3 2C PLEASANT A Bss1 Bss2 BC 2C1 2C2 NESS Ap BA Bw Bwss1 Bwss2 Bwss3 Bssk1 Bssk2 Bk BCk FETERITA Ap Bt1 Bt2 Bt3 Bt4 BC C LODGEPOLE 175 cm 150 cm 125 cm 100 cm 75 cm 50 cm 25 cm 0 cm

Problematic Clustering

# use local data
data("OSDexamples")

o <- reconcileOSDGeomorph(OSDexamples, 'hillpos')

knitr::kable(o$geom, digits = 3)
series Toeslope Footslope Backslope Shoulder Summit n shannon_entropy
AMADOR 0.000 0.000 0.667 0.111 0.222 36 1.224
ARGONAUT 0.062 0.125 0.458 0.312 0.042 48 1.856
CECIL 0.000 0.000 0.307 0.355 0.338 1027 1.582
DRUMMER 0.986 0.012 0.000 0.000 0.002 1617 0.116
HANFORD 0.517 0.329 0.154 0.000 0.000 149 1.436
MOGLIA 0.000 0.195 0.268 0.293 0.244 41 1.985
MUSICK 0.000 0.000 0.753 0.206 0.041 73 0.966
PALAU 0.200 0.200 0.200 0.200 0.200 100 2.322
PARDEE 0.058 0.000 0.192 0.212 0.538 52 1.650
PENTZ 0.013 0.000 0.625 0.112 0.250 80 1.357
REDDING 0.176 0.010 0.157 0.029 0.627 102 1.498
SIERRA 0.000 0.018 0.802 0.108 0.072 111 0.980
SYCAMORE 1.000 0.000 0.000 0.000 0.000 17 0.000
VLECK 0.000 1.000 0.000 0.000 0.000 2 0.000
WILLOWS 1.000 0.000 0.000 0.000 0.000 47 0.000
YOLO 0.593 0.322 0.085 0.000 0.000 59 1.275
ZOOK 0.986 0.015 0.000 0.000 0.000 482 0.109 
res <- vizHillslopePosition(o$geom, verbose = TRUE)
print(res$fig)

Proportion PARDEE REDDING CECIL MOGLIA PALAU PENTZ AMADOR MUSICK SIERRA ARGONAUT HANFORD YOLO DRUMMER ZOOK WILLOWS SYCAMORE VLECK 0.0 0.2 0.4 0.6 0.8 1.0 52 102 1027 41 100 80 36 73 111 48 149 59 1617 482 47 17 2 1.65 1.5 1.58 1.98 2.32 1.36 1.22 0.97 0.98 1.86 1.44 1.28 0.12 0.11 0 0 0 N H Toeslope Footslope Backslope Shoulder Summit 

# truncate profiles at 180cm
options(.aqp.plotSPC.args = list(max.depth = 180))

par(mar = c(1, 0, 1, 2))
plotGeomorphCrossSection(OSDexamples, type = 'line', maxIter = 100, j.amount = 0.05, verbose = TRUE)

A1 A2 Bt1 Bt2 2Bt3 2R PARDEE A1 A2 2Bt 3Bqm REDDING Ap Bt1 Bt2 BC C CECIL A Bkz Bknz 2Bknz 2Bknyz 3Bknyz MOGLIA A1 B1 B21 B22 B23 C PALAU A1 A2 Bw Bt Cr PENTZ A Bw1 Bw2 Cr AMADOR Oi A AB Bt1 Bt2 BCt MUSICK Ap BAt Bt1 Bt2 BCt C SIERRA A1 A2 Bt1 Bt2 Bt3 Cr ARGONAUT A1 C1 C2 HANFORD Ap1 Ap2 A1 A2 C1 C2 Ab C'3 YOLO Ap A BA Bg Btg1 Btg2 2Btg3 2Cg DRUMMER Ap A1 A2 A3 A4 Bg BCg Cg ZOOK Apg Bssg Bknssg Bnssg1 Bnssg2 Bnssg3 Bnssg4 WILLOWS Ap Bg1 Bg2 C SYCAMORE A1 A2 2Bt1/E 2Bt2 2Bt3/2Bq1 2Bq2 2Bq3 2Cr VLECK 175 cm 150 cm 125 cm 100 cm 75 cm 50 cm 25 cm 0 cm 0.0 0.2 0.4 0.6 0.8 1.0 Summit Shoulder Backslope Footslope Toeslope Proportion 

par(mar = c(1, 0, 1, 2))
plotGeomorphCrossSection(OSDexamples, type = 'line', clust = FALSE)

A1 A2 Bt1 Bt2 2Bt3 2R PARDEE Ap Bt1 Bt2 BC C CECIL A1 A2 2Bt 3Bqm REDDING A1 A2 Bw Bt Cr PENTZ A Bkz Bknz 2Bknz 2Bknyz 3Bknyz MOGLIA A Bw1 Bw2 Cr AMADOR Oi A AB Bt1 Bt2 BCt MUSICK Ap BAt Bt1 Bt2 BCt C SIERRA A1 A2 Bt1 Bt2 Bt3 Cr ARGONAUT A1 B1 B21 B22 B23 C PALAU A1 A2 2Bt1/E 2Bt2 2Bt3/2Bq1 2Bq2 2Bq3 2Cr VLECK A1 C1 C2 HANFORD Ap1 Ap2 A1 A2 C1 C2 Ab C'3 YOLO Ap A BA Bg Btg1 Btg2 2Btg3 2Cg DRUMMER Ap A1 A2 A3 A4 Bg BCg Cg ZOOK Ap Bg1 Bg2 C SYCAMORE Apg Bssg Bknssg Bnssg1 Bnssg2 Bnssg3 Bnssg4 WILLOWS 175 cm 150 cm 125 cm 100 cm 75 cm 50 cm 25 cm 0 cm 0.0 0.2 0.4 0.6 0.8 1.0 Summit Shoulder Backslope Footslope Toeslope Proportion 

## TODO: ensure geomorph data + type aren't mixed up

hydOrder(osd$hillpos, g = 'hillpos', clust = FALSE)
##  [1] "SPEARVILLE"   "HARNEY"       "RICHFIELD"    "ULYSSES"      "MANSKER"      "ULY"         
##  [7] "COLY"         "PENDEN"       "BUFFALO PARK" "DUROC"        "NESS"         "LODGEPOLE"   
## [13] "FETERITA"     "PLEASANT"
hydOrder(osd$hillpos, g = 'hillpos', clust = TRUE)
## $clust
## 
## Call:
## diana(x = daisy(x.prop, type = list(numeric = 1:n.prop)))
## 
## Cluster method   : NA 
## Distance         : euclidean 
## Number of objects: 14 
## 
## 
## $hyd.order
##  [1] "SPEARVILLE"   "HARNEY"       "RICHFIELD"    "ULYSSES"      "MANSKER"      "ULY"         
##  [7] "COLY"         "PENDEN"       "BUFFALO PARK" "DUROC"        "NESS"         "LODGEPOLE"   
## [13] "FETERITA"     "PLEASANT"    
## 
## $clust.hyd.order
##  [1] "SPEARVILLE"   "HARNEY"       "RICHFIELD"    "DUROC"        "ULYSSES"      "MANSKER"     
##  [7] "ULY"          "COLY"         "PENDEN"       "BUFFALO PARK" "LODGEPOLE"    "NESS"        
## [13] "FETERITA"     "PLEASANT"    
## 
## $match.rate
## [1] 0.3571429
## 
## $obj
## [1] 452
hydOrder(osd$flats, g = 'flats', clust = FALSE)
## [1] "MANSKER"   "ULYSSES"   "RICHFIELD" "FETERITA"  "DUROC"     "LODGEPOLE" "NESS"

Reset options.

options(.aqp.plotSPC.args = NULL)

This document is based on aqp version 2.2, soilDB version 2.8.8, and sharpshootR version 2.3.3.