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 9 slots
##  $ competing       :'data.frame':    61 obs. of  3 variables:
##  $ geog_assoc_soils:'data.frame':    75 obs. of  2 variables:
##  $ geomcomp        :'data.frame':    12 obs. of  9 variables:
##  $ hillpos         :'data.frame':    12 obs. of  8 variables:
##  $ mtnpos          : logi FALSE
##  $ terrace         :'data.frame':    6 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':    31 obs. of  5 variables:
##  $ pmorigin        :'data.frame':    5 obs. of  5 variables:
##  $ mlra            :'data.frame':    99 obs. of  4 variables:
##  $ ecoclassid      :'data.frame':    59 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':    24 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] 0.7857143
h2$match.rate
## [1] 0.7857143
# inspect objective function evolution
tr <- h2$trace

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

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

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)

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

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

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

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

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)

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

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)

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

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)

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

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)

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)

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)

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

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)

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.800 0.033 0.167 30 0.852
ARGONAUT 0.095 0.119 0.381 0.357 0.048 42 1.958
CECIL 0.000 0.000 0.307 0.355 0.338 1027 1.582
DRUMMER 0.988 0.010 0.000 0.000 0.002 1542 0.103
HANFORD 0.556 0.302 0.142 0.000 0.000 162 1.393
MOGLIA 0.000 0.195 0.268 0.293 0.244 41 1.985
MUSICK 0.000 0.000 0.694 0.282 0.024 85 1.008
PALAU 0.200 0.200 0.200 0.200 0.200 100 2.322
PARDEE 0.062 0.000 0.104 0.229 0.604 48 1.516
PENTZ 0.013 0.000 0.718 0.038 0.231 78 1.093
REDDING 0.198 0.011 0.154 0.022 0.615 91 1.502
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 24 0.000
YOLO 0.611 0.278 0.093 0.000 0.018 54 1.372
ZOOK 0.986 0.015 0.000 0.000 0.000 484 0.109 
res <- vizHillslopePosition(o$geom, verbose = TRUE)
print(res$fig)

# 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)

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

## 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"         "FETERITA"
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: 12 
## 
## 
## $hyd.order
##  [1] "SPEARVILLE"   "HARNEY"       "RICHFIELD"    "ULYSSES"      "MANSKER"      "ULY"         
##  [7] "COLY"         "PENDEN"       "BUFFALO PARK" "DUROC"        "NESS"         "FETERITA"    
## 
## $clust.hyd.order
##  [1] "ULYSSES"      "MANSKER"      "ULY"          "COLY"         "PENDEN"       "BUFFALO PARK"
##  [7] "SPEARVILLE"   "HARNEY"       "RICHFIELD"    "DUROC"        "NESS"         "FETERITA"    
## 
## $match.rate
## [1] 0.25
## 
## $obj
## [1] 350
hydOrder(osd$flats, g = 'flats', clust = FALSE)
## [1] "MANSKER"   "ULYSSES"   "RICHFIELD" "FETERITA"  "DUROC"     "LOFTON"    "NESS"

Reset options.

options(.aqp.plotSPC.args = NULL)

This document is based on aqp version 2.0.3, soilDB version 2.8.1, and sharpshootR version 2.2.